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INVERSION OF THE SPHERICAL MEAN 
TRANSFORM WITH SOURCES ON A HYPERPLANE 

ALEKSEI BELTUKOV 



Abstract. The object of this study is an integral operator S 
which averages functions in the Euclidean upper half-space R™ 
over the half-spheres centered on the topological boundary <9R™. 
By generalizing Norton's approach to the inversion of arc means 
in the upper half-plane, we intertwine S with a convolution oper- 
ator V . The latter integrates functions in R™ over the translates 
of a paraboloid of revolution. Our main result is a set of inversion 
formulas for V and S derived using a combination of Fourier anal- 
ysis and classical Radon theory. These formulas appear to be new 
and are suitable for practical reconstructions. 

1. Introduction 

The transform of the title S averages functions in the Euclidean up- 
per half-space = W 1 x (n > 1) over the half-spheres centered 
on the topological boundary <9IR™ +1 . Courant and Hilbert's well known 
result in jl] implies that S is injective on the space of continuous func- 
tions. If one further restricts the domain of S to rapidly decreasing 
functions, as done in pQ, then S can be shown to have a continuous 
inverse S" 1 . The main purpose of this article is to present a novel for- 
mula for that inverse which compliments the ones in [13], [6], pQ, [T3j . 
|15j . and [11] in several important ways: 

(1) Our inversion of spherical means holds in all dimension^, has 
explicit form involving the adjoint operator S* @, and applies 
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^The inversion formulas in |14j and [15] are derived only in dimensions two and 
three, respectively. 

2 The inversion formula in [T] is stated in Fourier domain. In practice it is often 
preferable to have a back-projection type formula (that is, involving the adjoint) 
as it leads to more stable reconstruction algorithms. 
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without restrictions to all rapidly decreasing functions in the 
Schwartz space & (K+ +1 ) I 

(2) The validity of our inversion formula follows directly from the 
classical Radon inversion. We thus establish a new relation 
between the spherical mean transform S and the classical Radon 
transform 1Z which, together with [5] and [2], brings the count 
for such relations to three. 

(3) Our formula for S~ l may be expressed as a composition of a 
certain convolution operator with a pair of simple one- 
dimensional mappings. That compositional structure lends it- 
self to a robust Fourier-based algorithm suitable for tomographic 
reconstruct ions . 

Although our interest in spherical means is purely mathematical, we 
would like to mention a few practical applications before delving into 
the theory. The following example, taken from |6j, explains why S 
continues to draw attention almost fifty years after having been first 
considered by Courant and Hilbert; it also explains why we refer to the 
centers of the half-spheres as sources. 

Consider the following nonhomogeneous hyperbolic equation with a 
delta-source at (x s ,y s ,0) G 1R 3 and zero initial conditions: 



d 2 u 



V u + f(x, y,z)u = S(x- x s )5(y - y s ) 6(z) 6(t), 



du 

u(x, y, z, 0) = — (x, y, z, 0) = 0. 



Think of (x s ,y s ,0) as a variable point in the xy -plane and denote 
by g(x s ,y s ,t) the solution of Problem (11 .11) evaluated at the location 
of the source: {x = x s ,y = y s ,z = 0}. The inverse scattering problem 
associated with Equation (II. ip is to determine the coefficient f(x, y, z) 
from g(x s ,y s ,t). This problem and its close variants arise in sonar 
and radar imaging (see [12] and [3], respectively), ultrasound diagnos- 
tics jH], and seismic tomography [6j. As explained in [Bj, by making 
the Born approximation one can interpret the data g(x s ,y s ,t) as the 
average of f(x, y, z) over the sphere of radius t/2 centered at (x s , y s , 0). 
If / is compactly supported in the upper (or lower) half-space, which is 
frequently the case □, the averages over spheres become averages over 



3 In [14] the spherical means are assumed to be band-limited in one variable; 
as shown in [llj . the inversion formulas in [5J, pQ and [13) require strong analytic 
constraints to be imposed on the functional domains in addition to rapid decrease. 

4 For instance, / could describe the acoustic reflectivity of underwater or un- 
derground environment both of which can be thought of as being semi-infinite in 
extent. 



INVERSION OF SPHERICAL MEANS 



3 



half-spheres. In that case the solution of the linearized inverse scat- 
tering problem for Equation fll.ip is equivalent to the inversion of the 
spherical mean transform S. 

The inspiration for our results comes from Norton's work in [H] on 
inversion of S in In Section [2] we will show how Norton's approach 
can be generalized to an arbitrary number of dimensions. Specifically, 
we will define a pair of mappings which intertwine S with an operator V 
acting on functions in the full Euclidean space. The latter operator is of 
convolution type: it integrates functions in IR n+1 over the translates of a 
paraboloid. We will exploit the convolutional nature of V in Section 
where we find V~ x using distribution theory. Our main result, the 
formula for S , follows from the formula for V~ x and has the general 
form 

f = y(S*o)C y oS)(f), /(x, Z /)G^(R^ +1 ), 

where fC y is a one-dimensional operator defined in Section [31 The 
advantage of our inversion formula over a similar result in [13] is that 
our one-dimensional operator K, is a simple pull-back of the classical 
A-operator, whereas the operator K in the formula / = (S* o K oS)(f) 
derived in [13] is a much more complicated fractional power of the 
Laplacian. Furthermore, as we show in Section [61 the presence of the 
A-operator, common in Radon theory, is not coincidental: our inversion 
formula for the convolution V follows directly from the classical Radon 
inversion formula. In fact, we show in Theorem 16.41 that the spherical 
mean transform S is, in a sense, more naturally inverted by the classical 
Radon transform. 

The author is deeply indebted to Dr. Eric Todd Quinto of Tufts Uni- 
versity for detailed and helpful suggestions concerning the preparation 
of the manuscript. He would also like to thank his colleague Dr. Jialing 
Dai for stimulating mathematical discussions of Radon and spherical 
mean transforms. 



2. Preliminaries 

Henceforth it will be convenient to identify R n+1 = M n x M} and, 
similarly, eR"xR| = {(x,y) \ x E M. n , y > 0}. The spaces of 

smooth compactly supported function in M. n+1 and will be denoted 
by C£° (M n+1 ) and (M" +1 ), respectively; for Schwartz spaces we will 
use the symbols S? (R n+1 ) and & 

Let / G y [n > 1) be the unknown function to be de- 

termined from its spherical means: g = S(f). We parameterize the 
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spherical mean data as follows: 

2 , r f( x + u , \/y 2 - \ u \ 2 ) 

(2.1) g( Xj y) = —- y l ~ n / — S >L 

v ; v ; \s n \ y J yV - \u\ 2 

\u\<y 

In Equation (12.11) \u\ = \/u\ + . . . + u\ is the Euclidean norm in M. n , 
du = du\ . . . du n is the Euclidean volume element, and 



(2.2) \S r 



27T 2 



r(f) 



is the well-known expression for the area of the unit sphere S" 1 C M n+1 . 
On occasion we will refer to the hyperplane (M™ +1 = R n where the 
half-spheres are centered as the centerset. 

The somewhat cumbersome choice of parametrization in Equation (12.11) 
is intentional: in fact, it is key to our strategy. Following Norton in [14"] . 
let us introduce the mappings Af : / i— ► F and M. : g i— > G by 

(2.3) F(«,») = ^. 

2/<0, 



(2.4) G(x,j/) 



H^lyV^a;,^), y >0, 



y < 0. 

Together Equations (I2.ip . (I2.3P and (12.41) imply the relation 

(2.5) G(x,y)= / F(x + u,y- \u\ 2 ) du, 



which we adopt as the definition of operator V. Furthermore, since 
the lower-case functions (/, g) are uniquely determined by the upper- 
case functions (F, G) via the inverse mappings A/" -1 : F i— ► / and 
: G i— ► g given by 

(2.6) = yF(x,y 2 ) \ y>0 , 

(2.7) ff (x,y) = ^y l - n G{x,y 2 )\ y>Q . 

Thus the spherical mean transform (12.11) and the convolution equa- 
tion (12. 5p are equivalent relations. 
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The operational meaning of Equations ( 12. ip . ( 12.31) . ( 12.41) . and ( 12.51) 
can be conveniently illustrated by the following commutative diagram: 

fey (R™ +1 ) — geS(y (R+ +1 )) 



C2.8) n 



M 



Fey (R n+1 ) — G eV {y (R n+1 )) 

For our purposes it is not necessary to characterize the spaces 5 fj^ (R™ +1 )) 
and P (y (R n+1 )). However, it will be useful to discuss the mapping 
properties of the convolution operator V on the space of compactly 
supported functions (R™ +1 ); we do that in Section HI 

As can be seen from Diagram 12.81 the operators S and V are in- 
tertwined by the mappings Af and M.. Consequently, inverting the 
operator S is equivalent to inverting the convolution operator V: 

(2.9) S- 1 =N~ 1 oV- 1 oM. 

We conclude this section by defining the adjoints V*, S* and a few 
additional operators needed to state our results in Section [3j 
As follows from Equation ( 12.51) . the L 2 -adjoint V* is given by 



(2.10) V*(G)(x,y)= G(x + u,y+\u\ 2 )du. 

In defining S*, however, we are going to follow [1] and [llffi 

(2.11) S*(g)(x,y)= [ g(x + u,^y 2 + \u\ 2 ) du. 



Finally, for a rapidly decreasing function <p(x,y) e y(W n+1 ), we de- 
fine the Hilbert transform H y (4>)(x,y), acting on the y- variable, and a 
related operator A y , by 



+oo 

(2.12) HM(x,y) = -(P.V.) / 

J y -t 



— oo 

n-1 



(2.13) Ky — {) ( ^°^' - = 1 ' 3 ,5,..., 

In Equations ( 12. 12ft and (12.131) the symbol d y denotes partial differen- 
tiation and (P. V.) stands for Cauchy principal value. 



5 Norton does not use the adjoint spherical mean transform in [14]; nor do Xu 
and Wang in [15) . Fawcett defines an "approximate inversion" operation which is 
quite similar to (|2.11[) . The adjoint S* used by Nessibi, Trimeche and Rachdi in [IB] 
differs from ours by a constant. 
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3. Main results 

Our main result — the inversion of the spherical mean transform S — 
follows from Theorem 13.11 stated below. 

Theorem 3.1. The following inversion formula holds: 

(3.1) F = 7T- n {V* oA y oV){F), Fe.y(W l+1 ). 

Equation (13.11) will be derived in Section [5] using Fourier analysis and 
the theory of distributions. However, the final proof of the validity of 
Theorem 13.11 will be deferred until Section [6] where it will follow from 
the classical Radon inversion formula. Meanwhile, we state and prove 
our main result. 

Theorem 3.2. The following inversion formula holds: 

(3.2) f = v(S*oKyoS){f), fey{w + +1 ) ) 

where the one- dimensional operator K y (acting on the y -variable) is 
defined by 

(3.3) M0)(^) = ^A t (tV0 (v /t)) 

t=y 2 

Proof. Combining Equations (12.91) and (13. ip . we obtain 

/ = 7T-" (AT" 1 oroA^oMo S)(f). 
We now define the mapping £ as a simple modification of Equation (12.31) : 

f(x,y/y), y>0, 
y < 0. 

Notice that the inverse mapping £~ l is given by 

(3.5) f(x,y) = F(x,y 2 ) \ y>0 , 

which is analogous to Equation ( 12. 6p . 

From Equations (I2.10p . (12. lip , and (13 Ah follows the intertwining 
relation 

(3.6) yS* =N~ l oV* oC. 

To apply Equation (13.61) . we insert the composition £ o £~ 1 into Equa- 
tion (13. 2p as follows: 

(3.7) / = 7T- n (TV -1 o V* o £ o CT X o A y o M o S)(f). 

Although £ o £~ 1 does not act as an identity, the validity of (13. 2p is 
unaffected because the values of V*(F)(x,y) for y > depend only on 
the restriction of F(x,y) to the upper half-space. Using the intertwin- 
ing relation (13.61) . we replace the triple composition (Af^ 1 oV* o£) with 



(3.4) £(/)(*, y) 
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yS*. Finally, by explicitly rewriting the composition (C" 1 o A y o Ai), 
using Equations (13.41) and (12.41) . we are led to the operator K, y defined 
by Equation (13. 3p . □ 

Given the practical importance of spherical means in low dimen- 
sions, we conclude this section with explicit inversion formulas for S in 
dimensions two and three. 

Corollary 3.3. Let g = S(f) where f is rapidly decreasing. The 
following inversion formula holds in M 2 ,: 



(3.8) /M= * fU.VofK^U 

K J-oc L JO 



y 2 + u 2 — s 2 



du, 



where the bracketed integral is interpreted in the sense of Cauchy prin- 
cipal value. 

In the analogue of Equation (13. 8j) can be written as 



f(x,y) 



y_ 

2tt 



g(x + u, yjy 2 + \u\ 2 ) 
{y 2 + \u\ 2 ) 3 * 



(3-9) 



g' y (x + u, y/y* 


+ 


U 


2 ) 


y 2 + 


u 


2 



9yy(x + u, ^Jy 2 + \u\ 2 ) 



du. 



u\ 



where now x = (xi, a^)- 



Proof. Differentiating Equation (12.41) (with n = 2) with respect to the 
^-variable, we obtain 



(d y G)(x,y) = G y (x,y) 



0, y < 0. 

Now the successive application of Hilbert transform (I2.12p leads to 



oc 

(H y G y )(x,y) = (P.V.) j 



g y (x, y/i) dt 



y-t 2Vi 



oc 

(P.V.)/ 



g y (x,s) 



ds. 



y 



Finally, applying the adjoint V* (see Equation (I2.10p ). we get the fol- 
lowing explicit relation between F and g 



F(x,y) 



1 

7T 



:i*.v. 



g y (x + v,s) 
y + v 2 — s 2 



ds 



dv, 
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from which the inversion formula (13.8p results after the substitution 

f(x,y) = yF(x,y 2 ) \ y>0 , 

in accordance with Equation (12.61) . 

The argument leading to Equation (13.91) is completely analogous and 
will therefore be omitted. The corollary is proved. □ 



4. Mapping properties of V 

In this section we examine the mapping properties of the convolu- 
tion operator V on the space of smooth compactly supported functions. 
The lemmas below do not fully characterize the range of V, yet they 
establish some useful facts that will guide us to the correct expres- 
sion for V~ l in Sections [5] and [6j We begin with a simple asymptotic 
estimate. 

Lemma 4.1. Let F(x,y) be a smooth compactly supported function 
in M™ +1 {n > 1) with support lying in the upper half-space . Fix 
yo) £ and consider 

G(v) = \V(F)(x + v,y + \v\ 2 )\, 

as a function of v £ IR n . There exist r > and h > (possibly depen- 
dent on (xo, yo) ) such that 

rnay p 

(4.1) lim \v\ G(v) < — \B n - l \ r"- 1 h, 

|?;|^oo ' 2 

where \B n \ denotes the volume of a unit ball in W 1 (\B°\ = 1). 

Proof. It is instructive to first prove the statement of the lemma for 
n = 1. 

Since F is compactly supported, we can choose r > and h > so 
that supp(F) is entirely contained within the rectangle [x — r, x + r] x 
[0,h\. 

As v goes to infinity, the parabola 

P = {(.x + v + u, y + v 2 - u 2 ) | u e R 1 }, 



passing through (x ,y ), will intersect only the horizontal sides of the 
rectangle at points A(^Jy + t> 2 , 0) and B(^y + v 2 — h, h), as shown 
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(x + v,y + v 2 ) 



(x - r, h) 



(x - r, 0) 




(x + r, 0) 



Figure 1. Estimation of G(v) in two dimensions 



in Figure [H Therefore 



G(v) 



F(x + v + u, yo + v 2 — u 2 ) du 



-y/yo+v 2 ~h 



F(xq + v + u, yo + v 2 — u 2 ) du 



< max |F| {^/vq + v 2 - y/y~ + v 2 - hj . 

It is now a simple matter to compute 

lim v G(v) = - max \F\ h. 

Notice that this is in agreement with Equation ( 14.1ft for n = 1. 

We now turn our attention to the case n > 1. Our strategy is going to 
be the same as in the preceding case, except that now the support of F 
will be confined to the interior of the cylinder C = B n (xo,r) x [0,h] 
with sufficiently large dimensions. In order to visualize the intersection 
of C with the paraboloid 



P = {Oo + v + u,y 



u e W 1 } 



we project both orthogonally onto M. n using the mapping ir : (x, y) G 



pn+l 



x G R n . The result is shown in Figure [2] below. 
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Figure 2. Estimation of G(v) in higher dimensions 



The lightly-shaded disk is the projection of the cylinder C. The large 
solid circle is the intersection of the paraboloid P with the plane y = 
containing the base of C . The dashed circle is the intersection of P 
with the plane y = h containing the top of the cylinder C. Of primary 
interest to us is the darkly-shaded quadrilateral bound by four circular 
arcs (representing spherical surfaces in M n ). This region is the orthog- 
onal projection tt(C fl P) which we now use to derive the following 
estimate: 



G{v) = l F{xq + v + u, y + \v\ 2 — \u\ 2 ) du 

/ F(xo + v + u,y + \v\ 2 — \u\ 2 ) du 

'tt(CT\P) 

< max|F| vol(7r(CnP)). 



To finish the proof we will find the volume of 7r(CnP). To this end, we 
express vol(7r(C HP)) as the difference of volumes of two lens-shaped 
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domains: 

vo1(tt(C n P)) = vol (B n {x , r) n B n (x + v, a/z/o + M 2 ) 

- vol (5 n 0o, r) n B n (x + v, y/y + \v\ 2 - h) 
'r 2 -y \ TT fr 2 -y + h 



(4.2) 



U 



2r \v\ 



+ (yo + \v\ 2 )%U 



(yo + H -h)*U 



U . 

2 r |u| 
2 |v | 2 + y — r 2 
2 M a/^/o + \v\ 2 / 

2 |f | 2 + y — r 2 — h 



2\v\ ^yo + \v\ 2 -h 



where 
(4.3) 



U{s) = \B n ~ x \ / (l-t 2 )— dt. 



The form of Equation (14.21) suggests the introduction of the following 
limits: 



L\ = r n lim \v\ 

\v\— too 



E7 



r 2 - y \ _ u (r 2 -y + h 
2r\v\ 



lim |v| (i/o + M ) 2 ^ 

|d|^oo 



2 r \v | 
2 |-u| 2 + y — r 2 
2 |v| y/yo + \v\ 2 



2 


V 


2 + yo- 


r 


2 -h 


2 


v\\/yo + 


V 


2 - h 



L 3 = lim \v | (y + \v\ 2 - h)* U 

\v\— >oo 

The required limit on the left-hand side of Equation (14.11) is now a 
simple linear combination of L\, L 2 and L 3 : 

lim \v\ G(v) = max \F\ {L t + L 2 - L 3 ). 

\v\—>oo 

To find Li, substitute \v I -1 = e. This leads to 



U 



Li = r" lim 



2r 



e)-U 



r' 2 -yo+h 
2r 



1 



e 2 

since Ef(0) = -Ifi^ 1 ]. 

Expanding U(s), given by Equation (14.31) . into a generalized series 
at s = 1, we obtain 



2^ 



= IB- 1 !— ^ (1-5)^ + ((i- fl 



n+3 
I 2 



12 



ALEKSEI BELTUKOV 



which, together with 



2 


\v 


\ 2 + yo- 




2 


V 


y/yo + 


v\ 2 



l-l-r 2 \v\- 2 + 0(\v\ 



leads to the values of Lo and L 3 : 



J 2 — -^3 

We conclude that, for n > 1 



2^ 

U = U = \B n - l \ r n+1 . 



n+1 



IT! r 

lim \v\ G(v) = max \F\ L x = !— L r "-i / lj 

|d|^oo 2 

as required. □ 

Remark 4.2. From Lemma [4.11 follows that the composition ("P* o P) 
cannot be freely applied to functions in (M n+1 ). Indeed, let F 6 
(M n+1 ) be nonnegative and let (xo,yo) e SU PP(-^)- Then 

{V*oV){F){x ,y )= [ V{F){x + v,y + \v\ 2 )dv^oo, 

since V{F)(xq + v,yo + \v\ 2 ) = 0(|f| _1 ). Accordingly, any inversion 
formula of the type F — (X oV* oV)(F) requires strong restrictions on 
the functional domain of V to assure the convergence of (V* a V)(F). 

The above analysis applies, with minor modifications, to the spher- 
ical mean transform S. As Klein shows in the inversion formulas 
of the type / = (X o S* o S)(f), derived in [6] and [T3J, can be in- 
terpreted conventionally only on a limited subspace of C c °° . To 
avoid convergence issues associated with the backprojection operation 
S* o S, Klein suggests to precede the adjoint S* with a derivative d y . 
Our solution, on the other hand, is to seek an inversion formula of the 
kind / = (S*oKoS)(f). 

Our next lemma characterizes the shape of the support of G — V(F); 
It will be particularly useful in Section 

Lemma 4.3. Let G = V(F) for some F £ C c °° (R n+1 ) (n > 1) and 
let II be any hyperplane in IR n+1 that meets the support of G in general 
position. The intersection supp(G) D II is compact. 

Proof. Figure [3] illustrates the lemma in two dimensions (n = 1). 

The dark square, say, [—a, a] 2 for some large a > 0, bounds the 
support of F: supp(F) C [—a, a] 2 . The lightly-shaded [/-shaped region 
is the infinite union of parabolas covering supp(G): 

supp(G) C (J {(x + v,y + v 2 ) | v G R 1 }. 
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supp(F) 



Figure 3. Intersection supp(G) n II in two dimensions 

Evidently, the intersection of the ^/-shaped region with a generic line II 
is compact. In fact, in two dimensions, any line meeting supp(G) forms 
a compact intersection. 

We now present a rigorous proof of the general case which we base 
on Heine-Borel Theorem. 

Since both II and supp(G) are relatively closed in M n+1 , the inter- 
section supp(G) nil is closed. Our task therefore is reduced to showing 
that supp(G) fl II is bounded. 

Let 7r : (x, y) G R" +1 i— ► x G M n denote the orthogonal projection, 
as in Lemma 14.11 Since the hyperplane II is in general position, the 
boundedness of supp(G) nil implies and is implied by the boundedness 
of the projection 7r(supp(G) fl II). We therefore direct our attention to 
the latter. 

To visualize 7r(supp(G) nil) we describe II explicitly with a standard 
equation u ■ x + 9y = p. Our assumption about II being in general 
position allows us to disregard the special cases u> = and 9 = 0. 
Hence we can normalize uo e IR n to have unit length and choose 9 G M 1 
to be positive. 

Now fix (x ,y ) G supp(F) and consider the intersection of II with 
the paraboloid of revolution: 

P = {{x + v,y + \v\ 2 ) | v G R n }. 

The projection ir(P fl II) is readily seen to be a sphere with radius 

= ^/l + 4:9(p-u-x -9y ) 
r 29 
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and center x + u/{2 9) G W 1 . Set 



W = 7TB W 1 + 46 * max 

2fc> Y (x,s/)esupp(F) 

We have the inclusion: 

UJ 



tt(p nn)c 5 n (x + 

Since supp(G) is covered by paraboloids 



2 ^ ) r max 



supp(G) C (J {(x + v,y + \v\ 2 ) | f G IR n }, 

(x-,y)esupp(F) 

the projection 7r(supp(Cr) PI II) is covered by balls: 

7r(supp(G) n n) C |J B n (x + r max ) . 

xe7r(supp(F)) 

We conclude that, since 7r(supp(F)) is compact and r max depends only 
on (no, 6), the projection 7r(supp(G)nlI) and, consequently, supp(G)nll 
must be bounded. This proves the lemma. □ 

Remark 4.4. Figure [3] shows that the function G(x, y) = V(F)(x, y) (in 
two dimensions) is not bandlimited in the y-variable, as assumed by 
Norton in [H]. Indeed, as a function of the y- variable alone, G(x,y) is 
compactly supported. Yet, as shown on Page 360 in [7], a compactly 
supported function cannot be bandlimited for its Fourier transform is 
analytic and cannot vanish outside a compact set. Accordingly, in our 
derivation of the inversion formula for V, we avoided any assumptions 
about the spectrum of G(x, y). As a result, our formula for V~ x is quite 
different from Norton's. 



5. Fourier inversion of the convolution V 

Our aim in this section is to derive Equation fl3.ll) which gives the in- 
version formula for the convolution operator V; recall that the latter is 
defined by Equation (12.51) in Section [2j The derivation largely amounts 
to straightforward de-convolution of Equation (12.51) which requires us 
to introduce Fourier notation and collect a few calculus facts. 

Let <j)(x,y) be an integrable function in M n+1 . We define the Fourier 
transform (/>(£, rj) by 

(5-1) $(£,v)= / exp(-i£-x-ir]y)4>(x,y)dydx. 
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Remark 5.1. The "hat" -symbol is often used to denote the classical 
Radon transform which will be introduced in Section To avoid any 
ambiguity, we henceforth agree to use the "hat" only as a shorthand 
for the classical Fourier transform in IR n+1 . 

To avoid explicit use of the Fourier inversion formula, we will need to 
recognize Fourier transforms of the operators introduced in Section [2j 
these we now present below. 

As follows from Equation (33.4) on Page 313 in (7J, with our choice 
of the constant in front of the integral in Equation (12.121) the Fourier 
symbol of the Hilbert transform 7i y is given by 



(5.2) H y = -i sgn(??) 



-i 1] > 0, 
+i 1] < 0, 



that is, 7i y (<f)) = TC y <p = — i sgn(r/) <p. 

Proposition (17.2.1) on Page 157 in [7j show!] that the Fourier sym- 
bol of the partial derivative = d y is given by d y — ir\\ (Of course, 
this simple result can be easily derived from scratch using integration 
by parts.) Together with Equation (15.21) . the formula for d y leads to 
the Fourier symbol of the A-operator: 

(5.3) K=\v\ n , n = 1,2,3,... 

The latter is a familiar fact in signal processing and tomography. 

It now remains to establish some calculus facts which we prefer to 
keep separate from Fourier manipulations. We shall require the follow- 
ing integrals found in [S] (Page 490 Equations (5) and (6)): 

r(fi) sin ^/i arctan | 
(/3 2 + «5 2 )t 

r(/i) cos (ji arctan | 

U 

Integrals (15 .4p and (15. 5p converge for 3?(/i) > 0, 3?(/9) > \^(S)\ and are 
needed in the proof of Lemma 15.21 below. 



oo 



(5.4) / x M exp(— (3 x) sin(5x) dx 



o 

oo 



(5.5) / x^ 1 exp(— (3x) cos(Sx)dx 



6 Gasquet and Witomski use the symmetric definition of the Fourier transform 
with a factor of 2ir in the exponential 
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Lemma 5.2. For any (el" and r] ^ 0, 

lim / exp (i £ • x — (e + irf) |x| 2 ) fix 

(5-6) "° , , 

7r \ 2 / |(| — -n n\r] 



exp i 
^1/ V 4 77 

Proof. To find the limit in Equation ( 15. 6ft . we first compute the integral 
using the following steps. 

Completing the square, we obtain 



exp I — ^ . ) / exp I — (e + ^77) 



4(e + i??) 
which upon substitution 



2(e + ir]) 



dx. 



x 



2(e + ir]) 
results in 



exp 



-(e-Hi,)l«l a d 



4(e + zr])_ 

We now switch to polar coordinates and integrate over the angular 
variable thereby obtaining 

IS"- 1 ! exp ( ^ f ex P (-(e + i7])r 2 ) r n ~ l dr, 

where IS 1 ™ -1 ) is the constant given by Equation (12. 2p . Now substitute 
r = y/s and apply Euler's formula exp(i 6) = cos(#) +i sin(0) to rewrite 
the integral as 

i^ n_i i ( \i? \ ( r < \ , 

exp / exp(— es) cos(ri s) s 2 ds 



J exp(— es) sin(i] s) s 2 ds^j 



Notice that the complex integral in Equation ( 15.61) is now reduced to 
two real-valued integrals similar to the ones listed in Equations (15.41) 
and (I5.5p . Setting (3 = e, 6 = r\ and \i = n/2, in tabular integrals (15.41) 
and (15.51) . and substituting the results into the preceding expression, 
we are led to 
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Now the limit as e — > + can be easily taken and the subsequent re- 
placement of the constant IS 1 ™ -1 ) with its value given by Equation (12.21) 
leads to Equation (15. 6p . as required. □ 

We are now ready to derive the inversion formula ( 13.11) . 

Derivation of Equation (13. ip . Applying the Fourier transform (15. ip to 
both sides of Equation (I2.5P and interchanging the order of integrations, 
we get the product relation G = V F in which the Fourier multiplier 



V 



\ exp(i £ ■ u — i 7] \u\ 2 ) du 

JR n 



is divergent for all n > 2. This necessitates a full switch to distributional 
perspective. 

Motivated by the form of Equation (12.51) . we define a linear func- 
tional T on the Schwartz space of rapidly decreasing functions 5? (M n+1 ) 
by 

^(0) = / ( x ! — \ x \ 2 ) dx. 

Evidently, T is linear and continuous on 5? (R n+1 ): as such, it is a 
tempered distribution in 5?' (M n+1 ) with a well-defined Fourier symbol 
fey (W n+1 ) . 

Next we reinterpret the right-hand side of Equation (12. 5p as the 
convolution of T e S? with a test-function F G (lR n+1 ) , i.e.: 

G = T * F. The result of application of the Fourier transform (15.11) to 
Equation (12.51) can now be restated as G = T F. Ostensibly, we have 
only replaced the symbol V with a new symbol T. Yet this subtle 
notational change makes all the difference. Indeed, the Fourier symbol 
V, having lost classical meaning, can now be identified with T — a well- 
defined tempered distribution in 5?' 

Following Kanwal in [10] , we use an exponential convergence factor 
to regularize T, i.e., we define 

T = lim / exp (i £ • x — (e + i rj) \x\ 2 ) dx, 

which is precisely the limit computed in Lemma [5. 2 iR 

7 The use of convergence factors had originated in mathematical physics as a 
heuristic way of computing "finite parts" of divergent integrals. Yet, it is now 
well-known that the "physics" regularization of divergent Fourier integrals is equiv- 
alent to the formal definition of the Fourier transform on S*' (R™) as long as the 
assumptions of the Lebesgue Dominated Convergence Theorem are satisfied. 

8 Lemma 15.21 (giving T) is largely identical to Kanwal's computation of the fun- 
damental solution of Shrodinger's equation in free space (Page 270 in [ID]). Kanwal 
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Using Equation ( 15. 6ft we conclude that 

whence follows the reciprocal relation F = (1/T) G which we write in 
the form 

It remains to identify the terms multiplying G on the right-hand side. 

The bracketed term, being the complex conjugate T corresponds to the 
adjoint T* (which in turn corresponds to V*) while the term |?7| n is the 
Fourier symbol of the operator A y given in Equation (15.31) . We are thus 
led to Equation (15.61) . as required. □ 

Remark 5.3. It may appear that along with Equation (13.11) we may 
form other inversion formulas for the operator V by permuting V* 
with the operators defining A y . However, our discussion of the map- 
ping properties of V in Section H] prohibits V* from being the leftmost 
operator (Remark 14. 2p . It is quite possible that the inversion formula 
F = n- n (A y oV*o V)(F) may be interpreted on & (R n+1 ) in distribu- 
tional sense. Yet, from the practical point of view, an expression with 
classical meaning is far more preferable. 

It remains to verify that the inversion formula given by Equation (13.11) 
is well-defined which will then complete the proof of Theorem 13.11 We 
do this in the next Section [6] by deriving Equation (13. ip from the clas- 
sical Radon inversion formula. 

6. Connection to Radon theory 

In this section we are going to adhere to our earlier convention 
whereby we identify M. n+1 = W 1 x K 1 . Accordingly, we shall param- 
eterize points on the unit sphere S n C R™ +1 using pairs (u, 9) with 
uj G M n and 9 e M 1 . This may, unfortunately, make some of the equa- 
tions listed below look unfamiliar. Nevertheless, these are all straight- 
forward adaptations of the formulas in classical Radon theory compiled 
from the first chapter of Helgason's text [H] . 



71 



71 



exp —i 



nn \r)\ 



4:7] 



further mentions that a faster way to compute the integral is to use the standard 
Fourier transform of a Gaussian. We, however, we prefer the method used in this 
section as more informative. 
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Let (j) E y (M n+1 ) be a rapidly decreasing function. We introduce 
the classical (n + l)-dimensional Radon operator 1Z by 



(6.1) K((j))((uj,6),p) = (j)(x,y)dm(x,y). 

J (uj,8)-(x,y)=p 

Notice that according to Equation ( 16. II) the Radon operator 7Z takes 
G (M n+1 ) into its integrals over the n-planes in R n+1 parameterized 
by (uj, 9) ■ (x, y) = p. In the parametrization of the planes the normal 
vector (a;, 9) G S n is the point on the unit sphere, split in accordance 
with our earlier explanation, while p G M 1 is the signed distance to the 
origin. The measure dm(x,y) in Equation (16.11) is the n-dimensional 
Euclidean volume element in the plane. 

Now let ip — 7^(0) ■ Notice that, according to Equation (16.11) . if) is a 
function on S" x I 1 . Following Helgason in [9], we define the adjoint 
Radon transform 7Z* by 



(6.2) TVW(x, y) = / rf>((u, 9), (u, 9) ■ (x, y)) d{u, 9), 

where d(u, 9) is the normalized Euclidean volume element on S n . 

In Section[2]we introduced the Hilbert transform H y and A^-operator 
acting on the y- variable (Equations (I2.12p and (12.131) . respectively). In 
this section we will use the same operators but acting on the p-variable 
instead: 

+oo 

1 m r <t>{->t) 



nM(;p) = ^(p-v.) J 



p — t 



-l)^(H p od;), n = 1,3,5, 
71 = 2,4,6, 



A p = 

As can be readily deduced from Chapter I of [9JJ3 the following com- 
mutation relation holds: 

(6.3) (K o A y )(0)((u;, 9),p) = 9" (A„ o n){cf>){{u, 9),p). 

Most importantly, A-operator occurs in the following variant of the 
classical Radon inversion formula 

(6.4) = c n (TV o A p o K)(<f)), c n = (-2 7r)- n ^ppy > 

which strongly resembles the inversion formula for the convolution op- 
erator V stated in Theorem 13.11 of Section [31 

9 Hel gason's definition of A-operator does not have the powers of (—1) while his 
definition of Hilbert transform, unlike ours, has an imaginary unit as a multiple. 
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We complete our review of classical Radon theory by taking note of 
the well-known translation property of the operator TZ. Let r denote a 
shift of a function as in r( a fe )(0)(x, y) = (p(x + a,y + b). Then by making 
a straightforward substitution in Equation (16.11) . it is easy to show (see 
Page 3 in [9]) that 

(6.5) (K o T ia>b) )(<f>)((Lu, 9),p) = K{f){(u, 6),p+ (u, 9) ■ (a, &)). 

Now, in order to deduce Theorem 13.11 from Equation (16.41) . which is 
the principle goal of this section, we are going to require the following 
technical lemma. 

Lemma 6.1. Let = {(u,8) G S n | 9 < 0} be the lower half of the 
unit sphere S n G M n+1 . Then 

(6-6) / hi ( !_ - d " =\S n \ [ h(u 7 8)d(u,8), 
Ju n \z\ z J (1 + \z\ z ) 2 Js™ 

where d(u,9) is the normalized spherical measure (restricted to S r l) 
while the constant \S n \ is the total n-dimensional volume of the unit 
sphere given by Equation (12.21) . 

The proof of Lemma ISTTl is a straightforward, albeit involved compu- 
tation; it is presented in Appendix |A] Meanwhile, we proceed to the 
proof of the main results of this section. 

Let Fey (W 1 ) be a rapidly decreasing function and set G = V(F) 
as in Equation (12. 5p . Theorem 16.21 below validates the inversion for- 
mulas stated in Theorem 13.11 of Section [3j 

Theorem 6.2. The inversion formula 

F = n- n (V* o A y )(G) = u- n (V* oA y o V)(F), 

stated in Theorem \3.1\ in Section^ is valid for rapidly decreasing func- 
tions Fey (R n ) . 

Proof. Since F is rapidly decreasing, the operators A y and V commute: 
(A y oV)(F)(x,y) = (VoA y )(F)(x,y)= [ A y (F)(x + u, y - \u\ 2 ) du. 



Hence, for a fixed (xo,yo) G W 1 x R 1 , we can rewrite the right-hand 
side of the inversion formula as 



71 



A y (F)(x + u + v, y + \v\ 2 — \u\ 2 ) du 



dv, 



which evidently is an integral in R n x 1". Now the change of variables 
v + u = w, v — u = z, (\v\ 2 — \u\ 2 — (v + u) • (v — u) = w ■ z) 
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transforms that integral into the expression 



(27T)-* 



A y (F)(x + w,y + z ■ w) dw 



dz. 



which we proceed to restate in terms of Radon transform (16.11) . To see 
how the operator 1Z enters the picture, rewrite the bracketed integral 
as 

1 



0(x ,yo) ° A y )(F)(w,z ■ w) a/1 + \z\ 2 dw. 



This shows that A y (F)(x,y) is first shifted by (x ,y ), then integrated 
over the plane (w, z • w,w G lR n ) with respect to the natural volume 
element + \z\ 2 dw, and, finally, scaled by (1 + |^| 2 ) -1 / 2 : in sum, 

(ft ° r ixom) o A y )(F) (± ( !_ , . 



v 1 + M V v 1 + \ z \ 

where we have complete freedom in choosing the sign of the normal — an 
important symmetry to be exploited later. 

Choosing the positive sign for now and using the translation prop- 
erty (16.51) combined with the A-commutation property (16.31) . we are led 
to 

(A , oR)(F) f_^) ,_^) . {xo , y o) y 



(1 + H 2 ) » ' ' ■ \V l + \*. 
which reduces our task to proving 

F( IOl!D ) = (-2,r \; ' V -<**■ 

JR" (1 + 2 

We claim that the latter is equivalent to the Radon inversion for- 
mula ( 16 Ah . To prove the claim and thus settle the theorem, we apply 
the mapping $ from Lemma 16.11 to change integration over M. n into 
integration over the lower half of the unit sphere SI : 

I On I r 

F( % t/o) = rVw / (A p o H)(F)((u, 9), (u, 9) ■ (x ,y )) d(u, 9). 
Js™ 

We now observe that the integration over the lower half-sphere S 1 ™ is 
the consequence of our earlier choice of the positive sign in front of the 
normal 

± (*,-!) 
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of the plane {w, z ■ w,w E M n } whereas had we chosen the negative 
sign, we would have obtained instead the same integral over the upper 
half-sphere S™ as immediately seen from Lemma ISTTl Therefore we can 
extend the domain of integration to the entire unit sphere S n : 

1 I S n I f 



2 (-2tt)« J s 
1 \S n 



2 -2vr 



(1VoA p oTl)(F)(xo,yo) 



Now substitution of the area of the unit sphere S n with its expres- 
sion (12. 2p from Section [2] leads to Radon inversion formula (I6.4p , as 
required. □ 

The intimate connection between V and 1Z elucidated by Theo- 
rem 16.21 suggests that instead of using Fourier transform to invert V, 
as we did in Section we could instead use Radon transform. More- 
over, the use of Radon transform is, in some sense, more natural^ 
Indeed, Lemma [4.31 in Section H] shows that, unlike Fourier transform, 
the Radon transform in K n+1 can be freely applied to G — V(F) with- 
out recourse to distribution theory. This motivates our derivation of 
a Radon-based inversion formula (I6.10p in Theorem 16.41 below. First, 
however, we need to introduce some additional background. 

Let G y x R 1 ) : think of as the Radon transform of a 
rapidly decreasing function in M n+1 . We now define a one- dimensional 
operator Q p acting on the p- variable by 

(6.7) Qp(0)((o;,0),p)= / <P((u;,6),p+(u,9)-(u,-\u\ 2 ))du. 
The adjoint operator Q* is then given by 

(6.8) Q;(<f>)((uj,e),p)= f 4>((u;,9),p+(u;,9)-(v,\v\ 2 ))dv, 

while the inverse Q~ l is presented in Lemma [6.31 below. 

Lemma 6.3. Let Q p be the operator defined by Equation (16.71) . Then 
\f)\ n 

(6.9) = LL( Q ;oA p )(0), ^y^xR 1 ), 
where Q* is the adjoint operator defined by Equation (16. 8p . 



10 See also the discussion of applications of Radon transform to partial differential 
equations at the end of Chapter I in [5]. 
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Proof. Applying one-dimensional Fourier transform in the p- variable 

oo 

0((w,0),p) ^ ?(0 = y exp(-i£p)0((w,0),p)dp, 

— oo 

to the convolution relation ip = Q p (<p), we get the product relation 
ip = Qp<f> where the Fourier multiplier Q p can be regularized with an 
exponential convergence factor as in Lemma 15.21 in Section [5j 

Q p = exp(i£ (oj, 9) ■ (u, — \u\ 2 )) du 

JM." 

= lim / exp(i^u ■ u — (e + i^9 \u\ 2 )) du. 
Now, as follows from Equation (I5.6p . 

n _ 



and the reciprocal relation <p = (1/ Q p ) ^ can be written as 
- \9\ n 



7T 



7T V / If ^1 - Trn If 01 
exp I — z " 



Recognizing the bracketed expression as Q p = Q* and |f | n as the fa- 
miliar Fourier symbol of the A-operator (acting on the variable p), we 
are led to the statement of the lemma. □ 

Using Lemma 16.31 and Radon inversion formula (16.41) . it is easy to 
invert V as Theorem 16.41 demonstrates. 

Theorem 6.4. Let F e ^(!R ri+1 ) . The following inversion formula 
holds 

(6.10) F = (R,- 1 oQ- 1 oKoV)(F), 

where Q^ 1 is given by Equation (16.91) in Lemma \6. 3[ 

Proof. Applying Radon transform (16. ip to Equation (12.51) . interchang- 
ing the order of integrations, and using the translational property (16.51) 
of Radon transform TZ, we obtain the relation 

K(G)((cu,9),p) = f K(F)((cu,9),p+(cu,9)-(u,-\u\ 2 ))du, 

which we can rewrite in operator form as 

K(G) = (Q p oK)(F). 
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Now writing G = V(F) and applying the inverses Q~ x and in 
that order, we get the statement of the theorem. □ 



7. Conclusions 

The inversion formulas for the spherical mean transform S, pre- 
sented in Sections [3] and strengthen the existing results and shed new 
light on the relation between spherical means and classical Radon the- 
ory. The Fourier-Radon techniques used to obtain our results have the 
virtue of simplicity, although this, admittedly, is a relative term. Unfor- 
tunately, the same techniques are unlikely to apply to spherical mean 
transforms with curved centersets. Indeed, from the group-theoretic 
point of view, our inversion of the operator S hinges on replacing its 
symmetry group with a much simpler symmetry group of the oper- 
ator V: notice that the symmetries are being replaced, not created. 
Yet, there are few centersets in IR n with sufficiently large symmetry 
groups. E3 In fact, flat centersets are the most symmetric, followed by 
spheres. The symmetry group of a spherical centerset SO(n) has lower 
dimension than that of V. This means that averages over the spheres 
centered on a sphere S 1 " -1 C W 1 cannot be straightforwardly related to 
any convolution in W 1 . Nevertheless, we consider our ad hoc inversion 
worthy of attention. The compositional structure of S" 1 exhibited in 
Equation (12.91) suggests a straightforward and effective FFT-based re- 
construction algorithm. We plan to pursue the development of such an 
algorithm and its numerical analysis elsewhere. 

It would also be of interest to extend the two Lemmas in Section H] 
to a full range characterization of spherical means S. We conjecture 
that a continuous function G which has asymptotic behavior 

G(x + v,y+ \v\ 2 ) = 0(M -1 ) 

given by Equation (14. ip is the image of a rapidly decreasing function F 
under V . If, as in Lemma I4.3[ the support of G can be covered by 
a union of parabolas with vertices ranging over a compact set, the 
pre-image F must be compactly supported. 



11 In R 2 all symmetric, or homogeneous, centersets can be enumerated by classi- 
fying the one-parameter subgroups of the Mobius group of the plane. The complete 
list, presented in the author's doctoral thesis, is limited to lines, rays, circles, open 
intervals, and logarithmic spirals. 
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Appendix A. Proof of Lemma 16.11 

Proof of Lemma \6.1\ For notational convenience, let us define mapping 
$ . R« s r l by 



Geometrically, $ projects the points of the tangent space T(o,...,o-i){S n ) 
onto as illustrated by Figure HI we may therefore call $ a central 
projection. 




Figure 4. The central projection $ 
Notice that $ is invertible with the inverse $ _1 : S 1 ™ i— > R n given by 

^((M)) = 

notice further that $ is smooth: it is therefore a diffeomorphism be- 
tween the Euclidean space R n and the unit half-sphere S™. 

We shall now introduce some additional terms from differential geom- 
etry. Let vols" denote the Riemannian volume form on the unit sphere 
S n C R n+1 (restricted to the lower half-sphere 5") and, likewise, let 
vo1r™(z) = dz± A . . . A z n denote the natural volume in R n ; we shall 
use the absolute value sign to distinguish a differential form from its 
corresponding density; finally, we denote by $* : A k (S™) i— > A fc (R n ) the 
pull-back mapping induced by $ : R n i— > S™ on the space of fc-forms. 
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We can now restate Equation (16.61) abstractly as 

VOl M n(z) 



(>o$)0) 



:i + \z\ 



n + l 



h((u,6)) \vo\ S n(u,6)\, 



noting that the proof of the lemma reduces to showing the equality of 
two densities 



|$*(vols»)(*)| 



VoW (z) 



:i+\z\ 



which will follow if (f>*(vol S n)(z) and (1 + \z\ 2 ) ( n+1 )/ 2 vol B «(z) are 
shown to agree up to sign. Equivalently, we must show that 

$*(vol 5 n)(» 1 



V0l R n(2) (1 

Therefore we define vol^™ explicitly by 



n + l 



vo\ S u(uj, 9) = Ui dw 1 A...du i ...A duo n A d# 



i=i 



+ 6 1 cL>i A ... A cLu n , 

where the empty dot over du>i indicates omission of the z-th term in 
the wedge product. Now from the definition of the pull-back of a 
mapping in differential geometry, we conclude that 



$*(volsn)(z) 

Voiron (z) 



n 

J2(-l) n+1+i coi(z) $* Uu x A... (Li 



A A dO 



+ 6(z) $* (dui A ... A du n ) = det(A), 
where A is the following (n + 1) x (n + 1) matrix 



dz\ ' ' ' ''••>■- ^ 



.4 



8z n 

89 

dz n 



We use the "empty dot" accent rather than the more customary "hat" from 
differential geometry to stay true to our promise in Section [5] to reserve the "hat"- 
symbol exclusively for the Fourier transform in M. n 
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with the entries given explicitly by 



dzj 
89 



1 



:i + ui 



1 + \z\ 2 - zf, i = j. 



as can be easily ascertained. The proof of the lemma is now reduced 
to a determinant computation: namely, we must show that 

det(A) = ■ 1 



(1 



n + l 



We will find the determinant of A in two steps. First, we introduce the 
matrix B by 



B 



1 + \z\ 2 - z 2 



-Z n Z\ 
Zl 



—Z\ z r , 



1 + \z\ 



-1 



noticing that the first n columns of A are the corresponding columns 
of B scaled by (1 + | | 2 ) 3/2 wm j e the last column in A is obtained 
from the last column of B via scaling by (1 + l^l 2 )" 1 ' 2 . Consequently, 
we have 

det ( A ) = V, , *3n±i det(S), 

(1 + \Z\ 2 ) 2 

and our task therefore is to compute det(B). To this end, we introduce 
a third and final matrix 



C 



Zl z n 



-Zl 



Z n Zl 
-Zl 



2 + \z\ 2 . 

whose connection to B is expressed by the relation 



B 



A + \z\ 2 ) l-C, 



where I is the (n + 1) x (n + 1) identity matrix. From the relation 
between B and C follows the relation between their determinants 

det(B) = -P c (l + \z\ 2 ), 

where Pc(l + \z\ 2 ) is the value of the characteristic polynomial of C 
at (1 + \z\ 2 ). Thus to find det(-B) we need the characteristic polynomial 
of C which prompts us to examine the eigenvalues of C . 



28 



ALEKSEI BELTUKOV 



It is straightforward to verify by direct computation that C has 
rank 2 with the following (n — 1) eigenvectors sharing the zero eigen- 
value: 



" -Z 2 ' 




" -23 " 




Zn 


Zl 















J 


Zl 


, . . . , 















Zl 


. 




o . 




o 



pn+1. 



while the two remaining eigenvectors of C 



Zl 

z 2 

Z 3 



■l±y/T+ 



pn+1 



have conjugate eigenvalues 1 + \z\ 2 =F y/l + \z\ 2 . We conclude that the 
characteristic polynomial of C is given explicitly by 

p c {\) = det(C — A I) = A"" 1 ((A - 1 - \z\ 2 ) 2 - 1 - \z\ 2 ) . 

Hence, 



det(B) = -p c (l + \z\ 2 ) = (l + \z\ 2 ) n , 



which implies 

det(A) 

as desired. 



2\n 



:i + \z\ 



to-(i + N') 



vl + Ul 



□ 
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